# ======================================================================
# Copyright CERFACS (February 2019)
# Contributor: Adrien Suau (adrien.suau@cerfacs.fr)
#
# This software is governed by the CeCILL-B license under French law and
# abiding  by the  rules of  distribution of free software. You can use,
# modify  and/or  redistribute  the  software  under  the  terms  of the
# CeCILL-B license as circulated by CEA, CNRS and INRIA at the following
# URL "http://www.cecill.info".
#
# As a counterpart to the access to  the source code and rights to copy,
# modify and  redistribute granted  by the  license, users  are provided
# only with a limited warranty and  the software's author, the holder of
# the economic rights,  and the  successive licensors  have only limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using, modifying and/or  developing or reproducing  the
# software by the user in light of its specific status of free software,
# that  may mean  that it  is complicated  to manipulate,  and that also
# therefore  means that  it is reserved for  developers and  experienced
# professionals having in-depth  computer knowledge. Users are therefore
# encouraged  to load and  test  the software's  suitability as  regards
# their  requirements  in  conditions  enabling  the  security  of their
# systems  and/or  data to be  ensured and,  more generally,  to use and
# operate it in the same conditions as regards security.
#
# The fact that you  are presently reading this  means that you have had
# knowledge of the CeCILL-B license and that you accept its terms.
# ======================================================================

"""Functions to compute the bounds on repetitions to ensure a given precision.

Two bounds are implemented:

  1. The analytic bound.
  2. The minimised bound.
"""

import math


def compute_analytic_bound(
    time: float,
    epsilon: float,
    hamiltonian_number_in_decomposition: int,
    max_spectral_norm_in_decomposition: float,
    order: int = 1,
) -> int:
    r"""Compute an upper-bound on :math:`r`, the number of time-steps to perform to \
    ensure the given error.

    This bound is the simplest one, but also the loosest one.

    :param time: The desired time of simulation.
    :param epsilon: The target precision.
    :param hamiltonian_number_in_decomposition: The number of hamiltonian matrices in
        the sum decomposition.
    :param max_spectral_norm_in_decomposition: An upper bound of the maximum spectral
        norm of the hermitian matrices in the sum decomposition.
    :param order: The order of the Trotter-Suzuki formula used.
    :return: an integer :math:`r` representing the number of time-slices that need be
        used during the hamiltonian simulation to guarantee the target precision
        analytically.
    """
    tau = (
        hamiltonian_number_in_decomposition
        * abs(time)
        * max_spectral_norm_in_decomposition
    )
    tmp = 2 * 5 ** (order - 1) * tau
    val = max(
        tmp,
        math.pow(tmp, (2 * order + 1) / (2 * order))
        * math.pow(math.e / (3 * epsilon), 1 / (2 * order)),
    )
    return max(int(math.ceil(val)), 1)


def compute_minimised_bound(
    time: float,
    epsilon: float,
    hamiltonian_number_in_decomposition: int,
    max_spectral_norm_in_decomposition: float,
    order: int = 1,
) -> int:
    r"""Compute an upper-bound on :math:`r`, the number of time-steps to perform to \
    ensure the given error.

    This bound is the a little bit more advanced than the analytical one but is still
    very loose.

    :param time: The desired time of simulation.
    :param epsilon: The target precision.
    :param hamiltonian_number_in_decomposition: The number of hamiltonian matrices in
        the sum decomposition.
    :param max_spectral_norm_in_decomposition: An upper bound of the maximum spectral
        norm of the hermitian matrices in the sum decomposition.
    :param order: The order of the Trotter-Suzuki formula used.
    :return: an integer :math:`r` representing the number of time-slices that need be
        used during the hamiltonian simulation to guarantee the target precision
        analytically.
    """

    tau = (
        hamiltonian_number_in_decomposition
        * abs(time)
        * max_spectral_norm_in_decomposition
    )
    factor_tau = 2 * 5 ** (order - 1) * tau

    def cost_func(
        repetitions: int, return_value_if_overflow: float = float("inf")
    ) -> float:
        x = factor_tau / repetitions
        if x > 1:
            return return_value_if_overflow
        else:
            return (repetitions / 3) * x ** (2 * order + 1) * math.exp(x)

    # First find lower & upper bounds for r
    r = 1
    while cost_func(r) > epsilon:
        r <<= 1  # r *= 2

    # If we are already inside the epsilon, return
    if r == 1:
        return 1

    # Then, binary search!
    low, up = (r >> 1), r
    while low != up:
        new = (low + up) // 2
        if cost_func(new) > epsilon:
            if low == new:
                return up
            low = new
        else:
            up = new
    return up


def compute_commutator_bound(
    time: float,
    epsilon: float,
    hamiltonian_number_in_decomposition: int,
    max_spectral_norm_in_decomposition: float,
    order: int = 1,
) -> int:
    r"""Compute an upper-bound on :math:`r`, the number of time-steps to perform to \
    ensure the given error.

    This bound is the most involved and the most efficient provable bound I know of.
    It has been introduced in :cite:`1711.10980v1`.

    :param time: The desired time of simulation.
    :param epsilon: The target precision.
    :param hamiltonian_number_in_decomposition: The number of hamiltonian matrices in
        the sum decomposition.
    :param max_spectral_norm_in_decomposition: An upper bound of the maximum spectral
        norm of the hermitian matrices in the sum decomposition.
    :param order: The order of the Trotter-Suzuki formula used.
    :return: an integer :math:`r` representing the number of time-slices that need be
        used during the hamiltonian simulation to guarantee the target precision
        analytically.
    """
    pass
